Week 02 - Introduction to Reading and Plotting Spatial Data
Introduction into GIS in R
There are a few key spatial packages available for Spatial Analysis in R, which you need to install as you progress through these exercises. The most basic are
sffor working with vector datarasterfor working with raster datargdalfor an R-friendly GDAL interface
The installation of the geospatial libraries GDAL, GEOS and PROJ.4 varies significantly based on operating system (and can be entailed). All of these are dependencies for sf, the R package that we will be using for spatial data operations throughout this course. Please follow the Option A: Local installation instructions for your operating system on this page under the heading GDAL, GEOS, and PROJ.4 and continue through the R Packages. Click on the linked pages to make your life easier.
Task 1: Reading vector data
The sf package, created by Edzer Pebesma and colleagues, has dramatically simplified reading vector spatial data into R.
In this exercise you will read in three shapefiles (parks, playgrounds, and forests one point file and two polygon files) and one geojson (shelters) using st_read(). If you’ve read in the files correctly, you will see a standard R data frame except it will show some header metadata about the file and you’ll see a special geometry column which we will discuss later.
Instructions
- Load the
sfpackage. - All your datasets reside in the ‘data’ folder.
- Import the
forestsshapefile (“forests.shp”). - Import the
playgroundsshapefile (“playgrounds4326.shp”). - Import the
parksshapefile (“parks.shp”). - Import the
sheltersgeojson (“shelters.json”). - Use the
head()function and identify the first few features of each layer.
# Load the sf package
___(___)
# Read in the forests shapefile
forests <- ___("../data/forests.shp")
# Read in the parks shapefile
parks <- ___(___)
# Read in the playgrounds shapefile
playgrounds <- ___(___)
# Read in the shelters json
shelters <- ___(___)
# View the first few features of all layers
___(___)Questions:
- How many features does each layer contain and what kind of geometry are they?
- What is the CRS value in these objects?
Solution
Reading layer `parks' from data source `C:\Users\Adela\Documents\RStudio\1_Teaching\cds-prep\data\parks.shp' using driver `ESRI Shapefile'
Simple feature collection with 47 features and 10 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 561480.2 ymin: 6210157 xmax: 581216 ymax: 6237475
projected CRS: ETRS89 / UTM zone 32N
Reading layer `forests' from data source `C:\Users\Adela\Documents\RStudio\1_Teaching\cds-prep\data\forests.shp' using driver `ESRI Shapefile'
Simple feature collection with 32 features and 15 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 566833 ymin: 6212338 xmax: 582568.3 ymax: 6236083
projected CRS: ETRS89 / UTM zone 32N
Reading layer `playgrounds4326' from data source `C:\Users\Adela\Documents\RStudio\1_Teaching\cds-prep\data\playgrounds4326.shp' using driver `ESRI Shapefile'
Simple feature collection with 51 features and 5 fields
geometry type: POINT
dimension: XY
bbox: xmin: 9.990331 ymin: 56.03292 xmax: 10.30607 ymax: 56.26783
CRS: NA
Reading layer `shelters' from data source `C:\Users\Adela\Documents\RStudio\1_Teaching\cds-prep\data\shelters.json' using driver `GeoJSON'
Simple feature collection with 19 features and 6 fields
geometry type: POINT
dimension: XY
bbox: xmin: 563459.7 ymin: 6212589 xmax: 578934.2 ymax: 6232898
projected CRS: ETRS89 / UTM zone 32N
Simple feature collection with 6 features and 6 fields
geometry type: POINT
dimension: XY
bbox: xmin: 570961.1 ymin: 6216613 xmax: 577431.9 ymax: 6232898
projected CRS: ETRS89 / UTM zone 32N
bookbar type navn skov
1 Ja Shelters Shelter i Gjellerup Skov <NA>
2 Ja Shelters Shelter i H<f8>rhaven i skoven <NA>
3 Ja Shelters Shelter i Lisbjerg gammel skov <NA>
4 Ja Shelters Shelter ved Moesg<e5>rd Strand <NA>
5 Ja Shelters Shelter i Mollerup Skov <NA>
6 Ja Shelters Shelter i H<f8>rhaven p<e5> bakken <NA>
mi_style mi_prinx
1 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 1
2 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 2
3 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 3
4 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 4
5 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 5
6 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 6
geometry
1 POINT (570961.1 6223407)
2 POINT (576307.6 6219521)
3 POINT (572490.7 6232898)
4 POINT (577431.9 6216613)
5 POINT (574609.9 6229691)
6 POINT (576455.1 6219324)
Simple feature collection with 6 features and 15 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 566833 ymin: 6223092 xmax: 576573 ymax: 6232903
projected CRS: ETRS89 / UTM zone 32N
bookbar type navn skovnr areal_m2 min_x min_y max_x
1 Ja Skove_store Åbyhøj Skov 290 32656 -224407 191537 -224289
2 Ja Skove_store Brabrand Skov 240 39231 -225585 191149 -225221
3 Ja Skove_store Årslev Skov 250 76629 -228610 191124 -228087
4 Ja Skove_store Gjellerup Skov 190 216967 -224931 190847 -224289
max_y oprettet_a oprettet_d rettet_af rettet_dat
1 192017 az25000 2019-05-07 12:02:50.657 az25000 2019-05-07 12:02:50.657
2 191359 az25000 2019-05-07 12:02:50.657 az25000 2019-05-07 12:02:50.657
3 191374 az25000 2019-05-07 12:02:50.657 az25000 2019-05-07 12:02:50.657
4 191485 az25000 2019-05-07 12:02:50.657 az25000 2019-05-07 12:02:50.657
mi_style mi_prinx
1 Pen (2, 1, 32768) Brush (2, 32768, 16777215) 1
2 Pen (2, 1, 32768) Brush (2, 32768, 16777215) 2
3 Pen (2, 1, 32768) Brush (2, 32768, 16777215) 3
4 Pen (2, 1, 32768) Brush (2, 32768, 16777215) 4
geometry
1 MULTIPOLYGON (((571023.4 62...
2 MULTIPOLYGON (((570135.1 62...
3 MULTIPOLYGON (((567335.9 62...
4 MULTIPOLYGON (((571154.8 62...
[ reached 'max' / getOption("max.print") -- omitted 2 rows ]
Simple feature collection with 6 features and 10 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 567350.2 ymin: 6210157 xmax: 581216 ymax: 6237475
projected CRS: ETRS89 / UTM zone 32N
bookbar type navn areal_m2 oprettet_a oprettet_d
1 Ja Parker Aarslev Møllepark 40022 az25000 2019-05-07 12:02:19.407
2 Ja Parker Bananen 5667 az25000 2019-05-07 12:02:19.407
3 Ja Parker Bavneparken 35611 az25000 2019-05-07 12:02:19.407
4 Ja Parker Bækvejsparken 19846 az25000 2019-05-07 12:02:19.407
5 Ja Parker Engdalgårdsparken 81476 az25000 2019-05-07 12:02:19.407
6 Ja Parker Forteparken 65751 az25000 2019-05-07 12:02:19.407
rettet_af rettet_dat
1 az25000 2019-05-07 12:02:19.407
2 az25000 2019-05-07 12:02:19.407
3 az25000 2019-05-07 12:02:19.407
4 az25000 2019-05-07 12:02:19.407
5 az25000 2019-05-07 12:02:19.407
6 az25000 2019-05-07 12:02:19.407
mi_style mi_prinx
1 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215) 1
2 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215) 2
3 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215) 3
4 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215) 4
5 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215) 5
6 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215) 6
geometry
1 MULTIPOLYGON (((567616.2 62...
2 MULTIPOLYGON (((581047.6 62...
3 MULTIPOLYGON (((567728.9 62...
4 MULTIPOLYGON (((574082.9 62...
5 MULTIPOLYGON (((575090.4 62...
6 MULTIPOLYGON (((575844.2 62...
Well done, now you should see how easy it can be to read in shapefiles and you got your first taste of what an sf object looks like.
Task 2: sf objects are data frames
Spatial objects in sf are just data frames with some special properties. This means that packages like dplyr can be used to manipulate sf objects. In this exercise, you will use the dplyr functions select() to select or drop variables, filter() to filter the data and mutate() to add or alter columns.
Instructions
- Load the
dplyrandsfpackages. - Use the
filter()function fromdplyron theparksobject to create a new data frame limited to parks greater than 5ha. - Use the
nrow()function on your new object to determine how many parks greater than 5ha are in the dataset. - Use the
select()function from dplyr to limit the variables in your over5ha dataset to justnavnandareal_m2and create a new data frame. - Use the
head()function to check which variables exist in your new data frame. Does the data frame only have thenavnandareal_m2columns (the answer is no, why)?
# Load the sf package
___
# ... and the dplyr package
___
# Use filter() to limit to over5ha parks
___
# Count the number of rows
___(over5ha)
# Limit to navn and areal_m2 variables
over5_lim <- over5ha %>% ___(navn, areal_m2)
# Use head() to look at the first few records
___(over5_lim)Solution
[1] 17
Great! You can see why the
sf package is so nice – your spatial objects are data frames that you can smoothly manipulate with dplyr. The number of parks over 5ha is 17 You may have noticed that when you used select the default is to keep the geometry column even if you didn’t explicitly list it as a column in select.
Task 3: Geometry is stored in list-columns
A major innovation in sf is that spatial objects are data frames. This is possible thanks, in part, to the list-column.
A list-column behaves, to a certain extent, like any other R column. The main difference is that instead of a standard value such as a single number, character or boolean value, each observation value in that column is a piece of an R list and this list can be as complex as needed. The list column allows you to store far more information in a single variable and sf takes advantage of this by storing all geographic information for each feature in the list.
In this exercise, you will convert the data frame to what’s called a tibble with tibble::as_tibble() (Note that dplyr::tbl_df() is now deprecated).
Instructions
- Load tidyverse in your workspace.
- Create a simple data frame
dfthat includes a single columnausingdata.frame(). - Add a list-column
bto your data frame with thelist()function. - Use
head()to look atdf. - Use
as_tibble()to convert the data frame to a tibble and print it to the console. This is just for cleaner printing. - Pull out the third observation from columns
aandbusingbaseR (you’ll need square brackets like[3]).
# Create a standard, non-spatial data frame with one column
df <- ___(a = 1:3)
# Add a list column to your data frame
df$b <- ___(1:4, 1:5, 1:10)
# Look at your data frame with head
___(df)
# Convert your data frame to a tibble and print on console
___(df)
# Pull out the third observation from both columns individually
df$___[___]
df$___[___]Solution
a b
1 1 1, 2, 3, 4
2 2 1, 2, 3, 4, 5
3 3 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
# A tibble: 3 x 2
a b
<int> <list>
1 1 <int [4]>
2 2 <int [5]>
3 3 <int [10]>
[1] 3
[[1]]
[1] 1 2 3 4 5 6 7 8 9 10
You now have a better sense of what a list column is. You can see how it can be used to store far more information in a single variable than other types of columns. These list-columns are how sf stores detailed geographic information on each feature in a single record. Converting the data frame to a tibble is not necessary but a tibble can provide a better print out of the object.
Task 4: Extract geometric information from your vector layers
There are several functions in sf that allow you to access geometric information like area from your vector features. For example, the functions st_area() and st_length() return the area and length of your features, respectively.
Note that the result of functions like st_area() and st_length() will not be a traditional vector. Instead the result has a class of units which means the vector result is accompanied by metadata describing the object’s units. As a result, code like this won’t quite work:
# This will not work
result <- st_area(parks)
result > 1e+05Instead you need to either remove the units with unclass():
# This will work
val <- 1e+05
which(unclass(result) > 1e+05)or you need to convert val’s class to units, for example:
# This will work
units(val) <- units(result)
which(result > val)
length(which(result > val))Instructions
- Check that your
sflibrary is still active and theforestsshapefile in memory - Compute the area of the forest units.
- Create a histogram of the areas using
hist()to quickly visualize the data spread. - Filter the forests object with
filter()and limit to forests withunclass(areas)> 200000 (areas greater than 20 hectares). - Can you plot the geometry and their names?
- Plot the geometry of the result with
plot()andst_geometry().
# Compute the areas of the forests
areas <- ___(forests)
# Create a quick histogram of the areas using hist
___(___)
# Filter to forests greater than 30000 (square meters)
big_forests <- ___ %>% ___(___(___) > 200000)
# Can you plot the forests with their names?
# Plot just the geometry of big_forests
___(___(big_forests))Solution
[1] "Gjellerup Skov" "Mollerup Skov"
[3] "Vestereng" "Tranbjerg Skov"
[5] "Vilhelmsborg Skov" "Skødstrup Skov"
[7] "Hørret Skov" "Brendstrup Skov"
[9] "Viby Høskov" "Fløjstrup Skov"
[11] "Kirkeskov" "Lisbjerg Skov, den gamle del"
[13] "Lisbjerg Skov, den nye del" "Thorsskov"
[15] "Hestehaven" "Skåde Skov"
[17] "Havreballe Skov" "Riis Skov"
[19] "Moesgård Skov" "Mariendal Skov"
Excellent! Computing geographic information for your vector layers can be done with functions like st_area() and st_length(). As you saw in this exercise, these functions produce a result that can be used in additional calculations but you need to be careful because the result is a units object that requires a little additional processing like using unclass().
Task 5: Plot vector spatial objects
The function for making a quick map/plot is a function you are already familiar with, plot(). You can, for example, type plot(my_data) to see your spatial object. The default, though, may not be what you want. The plot() function, when applied to sf objects, will create a set of maps, one for each attribute in your data. Instead, if you want to create a map of a single attribute you can extract that attribute using, as an example, plot(my_data["my_variable"]).
Frequently you just want to plot the raw geometry with no attribute color-coding (e.g., adding county boundaries to a map of points). For this, you can use the st_geometry() function to extract the geometry and plot the result. You can either create a new object or you can nest st_geometry() within the plot() function.
Often, you also want to plot multiple spatial objects together to see if and how they relate. To do this, use the plot()function twice in sequence, separated by ; with the second instance containing the add=TRUE argument plot( ,add=TRUE). You can use the col argument to differentiate each layer by color.
Instructions
- Use
plot()to plot theforestdata using all defaults. - Plot just the
areal_m2attribute of forests. - Create a new object that is just the geometry of the forests object with
st_geometry(). - Plot the geometry of the forests (the object you just created).
- Plot the geometry of the forests and the parks together, making the parks pink and the forests green.
Solution
Well done! Yes, these plots are not super pretty but you can’t beat plot() for a quick look using few keystrokes. And remember you can use plot(st_geometry(geo_object)) to plot just the geometry of your object.
Task 6: Reading in raster data
The term “raster” refers to gridded data that can include satellite imagery, aerial photographs (like ortophotos) and other types. In R, raster data can be handled using the raster package created by Robert J. Hijmans.
When working with raster data, one of the most important things to keep in mind is that the raw data can be what is known as “single-band” or “multi-band” and these are handled a little differently in R. Single-band rasters are the simplest, these have a single layer of raster values – a classic example would be an elevation raster where each cell value represents the elevation at that location.
Multi-band rasters will have more than one layer. An example is a color aerial photo in which there would be one band each representing red, green or blue light (RGB).
Instructions
- Load the
rasterpackage. - Load the elevation grid DNK_msk_alt.grd with the
rasterfunction and assign toelevationobject - Read in the orthophoto image with
brick()(this is multi-band raster called “Aarhus_1m.TIF”). - Use the
class()function to determine the class of each raster object you read in. - Use the
nlayers()function to determine how many bands/layers are in each object.
# Load the raster package
___(___)
# Read in the mound elevation single-band raster
elevation <- ____(___)
# Read in the orthophoto image multi-band raster
aarhus <- brick(_____)
# Get the class for the new objects
class(___)
class(___)
# Identify how many layers each object has
nlayers(___)
nlayers(___)Questions:
- What are the dimensions (number of columns and rows) of the
elevationraster? - What is the resolution of
aarhusorthophoto? How many layers does it contain and what do they represent?
Now you’ve learned how to read in single and multi-band rasters. You should have noticed, based on the nlayers() function, that the elevation object has a single layer and the aarhus object has four.
Solution
[1] "RasterLayer"
attr(,"package")
[1] "raster"
[1] "RasterBrick"
attr(,"package")
[1] "raster"
[1] 1
[1] 5
[1] 0.008333333 0.008333333
[1] 1 1
Task 7: Learn about your raster objects
Instead of storing raster objects in data frames, the raster package stores spatial data in specially designed R classes that contain slots where the data and metadata are stored. The data and metadata can be accessed using a suite of functions. For example, the spatial extent (the bounding box) of the object can be accessed with extent(), the coordinate reference system can be accessed with crs() and the number of grid cells can be determined with ncell().
Instructions
- You should have
rasterpackage loaded, and have theelevationlayer and ortophoto image layer for Aarhus in memory. - Use the
extent()function to get the extent of theelevationandaarhuslayer. - Use the
crs()function to get the coordinate reference system ofaarhusandelevation. - Use the
ncell()function to determine how many grid cells are in theelevationlayer and theaarhuslayer.
# Get the extent of the elevation and aarhus object
___(___)
# Get the CRS of the aarhus and elevation object
___(___)
# Determine the number of grid cells in both raster objects
___(aarhus)
___(elevation)Solution
class : Extent
xmin : 8
xmax : 15.3
ymin : 54.5
ymax : 57.8
class : Extent
xmin : 573334
xmax : 576668
ymin : 6223334
ymax : 6226668
CRS arguments:
+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs
CRS arguments: +proj=longlat +ellps=WGS84 +no_defs
[1] 11115556
[1] 346896
Great work! Although rasters are not stored as data frames, the metadata can easily be extracted using functions like extent(), crs() and ncell().
Question 3: It makes sense that the extents and number of cells for Danish national elevation data and Aarhus city would be different because their sizes diverge, but why are the units so different? What are these units?
Task 8: Plot your raster object
Similar to what you saw in the exercises related to vector objects it’s often useful to quickly look at a map of your raster objects with the plot() function.
The raster package has added useful methods for plotting both single and multi-band rasters. For single-band rasters or for a map of each layer in a multi-band raster you can simply use plot(). If you have a multi-band raster with layers for red, green and blue light you can use the plotRGB() function to plot the raster layers together as a single image.
Instructions
- Plot the
elevationraster with theplot()function, it is a single-band raster. - Plot the
aarhusobject with theplot()function, it is a multi-band raster. - Plot the
aarhusraster withplotRGB()to see all layers plotted together as a single image. Check out the optionalstretchargument in this function, and set it tolinearto brighten the image up a little - Try to plot the two images together (on top of one another). Can you do it?
# Plot the elevation raster (single raster)
___
# Plot the aarhus raster (as a single image for each layer)
___
# Plot the aarhus raster as an image
____Solution
Nice work! As you can see, the plot() function can be used to plot single layers while the plotRGB() function can be used to combine layers into a single image. Plotting two raster objects of different extents, resolution and CRS can not be done without additional wrangling
Task 9: Vector and raster coordinate systems
In order to perform any spatial analysis with more than one layer, your layers should share the same coordinate reference system (CRS) and the first step is determining what coordinate reference system your data has. To do this you can make use of the sf function st_crs() and the raster function crs().
When the geographic data you read in with sf already has a CRS defined both sf and raster will recognize and retain it. When the CRS is not defined you will need to define it yourself using either the EPSG number or the proj4string.
Instructions
- Ensure the packages
sfandrasterand the objectsforestsandplaygroundsandaarhusandelevationare loaded in your workspace. - Use
st_crs()to identify if a CRS exists and what it is for theplaygroundsandforestsobjects. - Use the
st_crs()function to define/assign a CRS to theplaygroundsobject, utilising the EPSG number 4326 as that is the original projection (now lost). - Use
crs()to identify if a CRS exists and what it is for theaarhusandelevationobjects. - Do not use the
crs()function to define/assign a CRS to theelevationobject as it already has a CRS of its own and renaming it won’t change the properties of the file.
Solution
# Determine the CRS for the elevation and playgrounds vector objects
st_crs(forests)Coordinate Reference System:
User input: ETRS89 / UTM zone 32N
wkt:
PROJCRS["ETRS89 / UTM zone 32N",
BASEGEOGCRS["ETRS89",
DATUM["European Terrestrial Reference System 1989",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4258]],
CONVERSION["UTM zone 32N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",9,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["unknown"],
AREA["Europe - 6°E to 12°E and ETRS89 by country"],
BBOX[38.76,6,83.92,12]],
ID["EPSG",25832]]
st_crs(playgrounds)Coordinate Reference System: NA
# Assign the CRS to playgrounds - its CRS is actually 4326, it just got lost in
# transport
crs_1 <- 4326
st_crs(playgrounds) <- crs_1
# Determine the CRS for the aarhus and elevation rasters
crs(aarhus)CRS arguments:
+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs
crs(elevation)CRS arguments: +proj=longlat +ellps=WGS84 +no_defs
# Assign the CRS to elevation
crs_2 <- "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"Nice! You can determine what the CRS is using st_crs() for vectors or crs() for rasters. If one doesn’t exist, you can also use those functions to define the CRS.
Task 10:Transform your layers to a common CRS
In the previous exercise, when you ran st_crs() and crs() you may have noticed that the CRS’ were different for the different layers. The raster layer’s CRS began with +proj=longlat and the vector layer’s began with +proj=utm. In order to use these layers together in spatial analysis we will want them to have the same CRS.
In this exercise you will transform (sometimes this is called “project”) the objects so they share a single CRS. It is generally best to perform spatial analysis with layers that have a projected CRS (and some functions require this). To determine if your object has a projected CRS you can look at the first part of the result from st_crs() or crs() – if it begins with +proj=longlat then your CRS is unprojected.
Note that you will use method = "ngb" in your call to projectRaster() to prevent distortion in the elevation image.
Instructions
- Use the
crs()function with argumentasText = TRUEon theaarhuslayer to get the CRS as a string and save this asthe_crs. - Use
st_transform()to transform the vectorplaygroundsobject to the CRS inthe_crs. - Use
projectRaster()to transform the rasterelevationobject to the CRS inthe_crs. This will take a few seconds. - Use
st_crs()onplaygroundsandelevationto confirm that they now have the same CRS. They should all have a CRS that starts with+proj=utm.*
Solution
# Get the CRS from the aarhus object
the_crs <- crs(aarhus, asText = TRUE)
# Project playgrounds to match the CRS of aarhus
playgrounds_crs <- st_transform(playgrounds, crs = the_crs)
# Project to match the CRS of aarhus
elevation_crs <- projectRaster(elevation, crs = the_crs, method = "ngb")
# Look at the CRS to see if they match
st_crs(playgrounds_crs)Coordinate Reference System:
User input: +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs
wkt:
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["Unknown based on GRS80 ellipsoid",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1],
ID["EPSG",7019]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 32N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",9,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16032]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
crs(elevation_crs)CRS arguments:
+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs
Great work! This may be the least fun part of spatial analysis but knowing how to do it will save you from a lot of frustration later.
Task 11: Plot vector and raster together
If the layers do not share a common CRS they may not align on a plot. To illustrate, in this exercise, you will initially create a plot with the plot() function and try to add two layers that do not share the same CRS. You will then transform one layer’s CRS to match the other and you will plot this with both the plot() function and functions from the tmap package.
Note that for this exercise we returned all the layers to their original CRS and did not retain the changes you made in the last exercise.
With the plot() function you can plot multiple layers on the same map by calling plot() multiple times. You’ll need to add the argument add = TRUE to all calls to plot() after the first one and you need to run the code for all layers at once rather than line-by-line.
Make sure sf, raster and tmap are loaded in your workspace.
Instructions
- Try to plot the
playgroundsobject on top of theaarhusobject withplotRGB(aarhus)followed byplot(playgrounds, add = TRUE). Do you see the playgrounds? If not, use thest_transformto project to shared CRS. - Re-run the
plot()code from the instruction above, color theplaygroundsgreen and make the points stand out larger by setting thelwdargument to 3. You can also increase the transparency ofaarhus, by setting itsalphaat around 100. - Run the given
tmapcode. Note thattm_rgb()is used for multi-layered raster
Solution
# Plot aarhus and playgrounds (run both lines together) Do you see the
# playgrounds?
plotRGB(aarhus, stretch = "lin", alpha = 100)
plot(playgrounds_crs, col = "green", lwd = 3, add = TRUE)# # See if aarhus and playgrounds share a CRS st_crs(playgrounds) crs(aarhus) #
# Save the CRS of the aarhus layer the_crs <- crs(aarhus, asText = TRUE) #
# Transform the playgrounds CRS to match aarhus playgrounds_crs <-
# st_transform(playgrounds, crs = the_crs) # Re-run plotting code (run both lines
# together) # Do the playgrounds show up now? plot(aarhus);plot(playgrounds_crs,
# add = TRUE)
# Simply run the tmap code
library(tmap)
tm_shape(aarhus) + tm_rgb() + tm_shape(parks) + tm_polygons(col = "green") + tm_shape(playgrounds_crs) +
tm_dots(col = "yellow", size = 1) Great work! As you noticed, you mostly can’t plot layers together if they don’t have the same CRS. You’ll see later that there are exceptions but it is definitely best practice to ensure the layers you’ll work with have a common, projected CRS.
Bonus: Loading data from the internet
There are a lot of online resources for spatial data, such as OSM, AirBnB, GoogleMaps, DivaGIS, and GADM database. In the exercise you will learn to work with rasters from the GADM database which contains spatial data - vector and raster - for all current countries and is well-integrated with R. You can load the data directly into R with the function getData() and start processing. The available data are:
- SRTM 90 (elevation data with 90m resolution between latitude -60 and 60)
- World Climate Data (Tmin, Tmax, Precip, BioClim)
- Global adm. boundaries (different levels)
In the case of GADM you must also provide the level of administrative subdivision (0=country, 1=first level subdivision)
Instructions
- To read elevation with the
getData()function from GADM, select the following arguments:nameisaltfor altitude (elevation),countryneeds a 3 letter ISO code for Denmark; getData(‘ISO3’) let’s you see the codes,maskneeds to be set to TRUE as we wish to mask surrounding countries.
- to read administrative data, choose the following:
nameisGADMfor global administrative data,countryneeds a 3 letter ISO code for Denmark; getData(‘ISO3’) let’s you see the codes,levelshould be set to the level of administrative subdivision (0=country border, 1=first level subdivision, eg. Midtjylland region, 2 = municipalities,… ).